In [44]:
from __future__ import division, print_function
# Stdlib
import os
# Third-party
from astropy.io import fits
import astropy.coordinates as coord
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# Custom
import gary.coordinates as gc
import gary.dynamics as gd
import gary.integrate as gi
import gary.io as io
import gary.potential as gp
from gary.units import galactic
First I need to find an orbit that passes near the Sun with a significant distance gradient
In [2]:
potential = gp.SphericalNFWPotential(v_c=0.22, r_s=8, units=galactic)
In [3]:
w0 = np.array([0,5,2,-0.35,-0.05,0.])
nparticles = 1000
sat_mass = 2.5e5 # solar masses
In [4]:
t,w = potential.integrate_orbit(w0, dt=-1., nsteps=3000)
In [5]:
fig = gd.plot_orbits(w, marker=None)
Now integrate a ball of orbits to make a pseudo-stream
In [6]:
E_scale = (sat_mass / potential.mass_enclosed(w0[np.newaxis, :3].copy()))**(1/3.) * 0.65
rtide = E_scale * np.sqrt(np.sum(w0[:3]**2))
vdisp = E_scale * np.sqrt(np.sum(w0[3:]**2))
sigma = np.array(list(rtide)*3 + list(vdisp)*3)
In [7]:
(rtide[0]*u.kpc).to(u.pc), (vdisp[0]*u.kpc/u.Myr).to(u.km/u.s)
Out[7]:
In [8]:
leading = w0.copy()
leading[:3] -= leading[:3]/np.linalg.norm(leading[:3]) * rtide
trailing = w0.copy()
trailing[:3] += trailing[:3]/np.linalg.norm(trailing[:3]) * rtide
In [9]:
ball_w0 = np.vstack((np.random.normal(leading, sigma/1, size=(nparticles,6)),
np.random.normal(trailing, sigma/1, size=(nparticles,6))))
In [10]:
t,ball_w = potential.integrate_orbit(ball_w0, dt=-1., nsteps=3400)
In [11]:
fig = gd.plot_orbits(ball_w[-1,:nparticles], marker='.', linestyle='none', color='r')
fig = gd.plot_orbits(ball_w[-1,nparticles:], marker='.', linestyle='none', color='b',
axes=fig.axes)
# fig.axes[0].set_xlim(-30,30)
# fig.axes[0].set_ylim(-30,30)
# fig.axes[0].set_xlim(-1,1)
# fig.axes[0].set_ylim(4,6)
In [12]:
obs = gc.gal_xyz_to_hel(ball_w[-1,:,:3].T.copy()*u.kpc)# .transform_to(coord.ICRS)
plt.plot(obs.b, obs.distance, linestyle='none')
Out[12]:
In [72]:
def make_data(noise_level):
path = os.path.join(os.path.abspath("../data/"), "noise{}".format(noise_level))
if not os.path.exists(path):
os.mkdir(path)
distance_bins = np.arange(5, 35, 3.)
lbins = np.linspace(0,150,150)
bbins = np.linspace(0,65,65)
pad = 25
size = (len(bbins) - 1 + 2*pad, len(lbins) - 1 + 2*pad)
fig,axes = plt.subplots(3,3,figsize=(15,15),sharex=True,sharey=True)
H,xedge,yedge = np.histogram2d(obs.l, obs.b, bins=(lbins,bbins))
scale = 1000 / H.max()
i = 0
for l,r in zip(distance_bins[:-1],distance_bins[1:]):
ix = (obs.distance.kpc > l) & (obs.distance.kpc < r)
H,xedge,yedge = np.histogram2d(obs.l[ix], obs.b[ix], bins=(lbins,bbins))
counts = np.random.poisson(noise_level, size=size)
counts[pad:-pad, pad:-pad] = counts[pad:-pad, pad:-pad] + H.T * scale
# axes.flat[i].pcolor(xedge, yedge, counts, cmap='Greys_r',
axes.flat[i].pcolor(counts, cmap='Greys_r',
vmin=noise_level-5*np.sqrt(noise_level),
vmax=noise_level+5*np.sqrt(noise_level))
# axes.flat[i].plot(obs.l[ix], obs.b[ix], linestyle='none', color='w')
axes.flat[i].set_title("{:.0f} < d < {:.0f}".format(l,r))
i += 1
hdu = fits.PrimaryHDU(counts)
hdu.writeto(os.path.join(path,"dist{:.0f}_{:.0f}.fits".format(l,r)), clobber=True)
axes.flat[0].set_xlim(0,150)
axes.flat[0].set_ylim(0,65)
fig.savefig(os.path.join(path, "dist_cuts.png"))
In [74]:
for nl in 2**np.arange(7,16+2,2):
make_data(nl)
In [ ]: